This Markdown does the following:
For the moment the data is stored locally. We will further put a link to download the data.
setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")
library(sfnetworks)
library(dodgr)
library(sf)
library(sp)
library(spdep)
library(geoR)
library(spatstat)
library(tidygraph)
library(data.table)
library(dplyr)
library(purrr)
library(TSP)
library(openxlsx)
library(raster)
library(stars)
library(terra)
library(ggplot2)
library(ggspatial)
library(knitr)
library(tidyr)
library(cartography)
viajes_dist <- read.xlsx("Tables produites/19-12/ViajesEODH_dist_corr.xlsx") # Trip table
UTAM <- st_read(dsn = "Data", layer = "EMU2019_area") %>% st_transform(4326) # UTAM boundaries
## Reading layer `EMU2019_area' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N
Manzana <- st_read(dsn = "Data", layer = "Manzana_Hogares") # Household manzana
## Reading layer `Manzana_Hogares' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 21828 features and 56 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 569999.6 ymin: 493719.5 xmax: 625763.1 ymax: 557019.5
## Projected CRS: WGS 84 / UTM zone 18N
EF <- read.xlsx("Data/facteurs d'émissions.xlsx") # emission factors
limites_bogota <- st_read(dsn = "Data", layer = "Localidad_Municipio_2017") %>% # Municipalities boundaries
st_transform(4326)
## Reading layer `Localidad_Municipio_2017' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 564193.1 ymin: 483995.6 xmax: 634113.2 ymax: 570281.8
## Projected CRS: WGS 1984 UTM, Zone 18 North, Meter
Transmi <- st_read(dsn = "Data", layer = "Transmilenio") # Transmilenio network
## Reading layer `Transmilenio' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 588057.1 ymin: 503694.6 xmax: 606100.3 ymax: 527398.5
## Projected CRS: WGS 84 / UTM zone 18N
centroides_limites_bogota <- st_centroid(limites_bogota) %>%
mutate(lon = st_coordinates(.)[,1],
lat = st_coordinates(.)[,2])
centroides_limites_bogota$MUNI <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$MUNI[centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"] <- ""
centroides_limites_bogota$LOCA <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$LOCA[centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"] <- ""
centroides_limites_bogota$DISPLAY <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$DISPLAY[!((centroides_limites_bogota$NOMBR_MUNI == "BOGOTA" & centroides_limites_bogota$NOMBRE %in% c("BOSA", "SANTA FE", "CIUDAD BOLIVAR", "TEUSAQUILLO", "FONTIBON", "ENGATIVA", "SUBA", "USAQUEN", "USME", "SAN CRISTOBAL", "CHAPINERO", "KENNEDY")) | (centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"))] <- ""
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SUBA"] <- 4.719
centroides_limites_bogota$lon[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- -74.027
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "USAQUEN"] <- 4.71
centroides_limites_bogota$DISPLAY[centroides_limites_bogota$NOMBRE == "CIUDAD BOLIVAR"] <- "CIUDAD \nBOLIVAR"
The same zoom level is chosen for Bogotá and Lima for comparison purpose
zoom_level <- 9.5
lon_span <- 360 / 2^zoom_level
lat_span <- 360 / 2^zoom_level
zoom_to_Bogota <- c(-74.12, 4.72) # Bogota
lon_bounds_Bogota <- c(zoom_to_Bogota[1] - lon_span / 2, zoom_to_Bogota[1] + lon_span / 2)
lat_bounds_Bogota <- c(zoom_to_Bogota[2] - lat_span / 2, zoom_to_Bogota[2] + lat_span / 2)
UTAM$Bogota <- 1
contour <- UTAM %>%
group_by(Bogota) %>%
summarise(geometry = st_union(geometry))
Dist_UTAM <- viajes_dist[,c(1,9,18)] %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
bks = c(100000,200000,500000,1000000,2500000)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30", labels = c("100000", "200000", "500000", "1000000", "2500000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Total daily PKT per \nUTAM of residence", fill = "Total passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
filter(p17_Id_motivo_viaje == 1) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_work = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
bks = c(100000,200000,500000,1000000)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_work))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30", labels = c("100000", "200000", "500000", "1000000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Commuting to work", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
filter(p17_Id_motivo_viaje == 3) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_school = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
bks = c(10000,20000,50000,100000,250000)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_school))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Going to school \nor university", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,6,9,18)] %>%
filter(p17_Id_motivo_viaje %in% c(12,16)) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_leisure = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
bks = c(5000,10000,20000,30000,60000)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_leisure))+ geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "#bb2b30") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Leisure or physical activity", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
filter(modo_principal == "A pie") %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_walk = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
#bks = c(0.5,0.8,1.2,1.5,2)
bks_walk <- round(as.numeric(quantile(UTAM$dist_walk/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_walk
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_walk/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Walking", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
filter(modo_principal == "Auto") %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_auto = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
UTAM$dist_auto[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA
#bks = c(1,2,3,4,5)
bks_car <- round(as.numeric(quantile(UTAM$dist_auto/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_car
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_auto/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82", na.value = "grey80") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Private car", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
filter(modo_principal %in% c("Transporte informal", "Bicitaxi")) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_informal = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
bks = c(0.2,0.5,1,2,3)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_informal/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82", na.value = "grey80") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Informal transport", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
filter(modo_principal %in% c("TransMilenio", "Alimentador", "Cable")) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_BRT = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
UTAM$dist_BRT[UTAM$UTAM == "UTAM2"] <- NA
bks = c(1,2,4,6,8)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_BRT/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_sf(data = Transmi, aes(color = "BRT (Transmilenio)"))+
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82",na.value = "grey80") +
scale_color_manual(values = "red")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Bus Rapid Transit (BRT)", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Dist_UTAM <- viajes_dist[,c(1,9,10,18)] %>%
filter(modo_principal %in% c("TransMilenio", "Alimentador", "Cable", "Intermunicipal", "SITP Provisional", "SITP Zonal")) %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(dist_transit = sum(f_exp*dist_real/1000))
UTAM <- UTAM %>% left_join(Dist_UTAM, by = "UTAM")
UTAM$dist_transit[UTAM$UTAM == "UTAM2"] <- NA
#bks = c(2,4,6,8,10)
bks_transit <- round(as.numeric(quantile(UTAM$dist_transit/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_transit
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = dist_transit/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#f8f8f8", high = "#103D82",na.value = "grey80")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Public transport", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM <- UTAM %>% st_transform(32618)
UTAMCentroids <- st_centroid(UTAM)
listPPV <- knearneigh(UTAMCentroids, k = 1, longlat = TRUE) # Finding each UTAM nearest neighbor
PPV <- knn2nb(listPPV, row.names = UTAM$UTAM) # Coercing knn objects into nb objects
distPPV <- nbdists(PPV, UTAMCentroids) # Computing the distance between nearest neighbors
print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 553.5302 1365.71 1587.616 2181.223 1854.208 10899.47
hist(unlist(distPPV), breaks = 20,
main = "Distance to the closest neighbor",
col = "black", border = "white", xlab = "Distance", ylab = "Frequence")
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
summary(nb)
## Neighbour list object:
## Number of regions: 132
## Number of nonzero links: 632
## Percentage nonzero weights: 3.627181
## Average number of links: 4.787879
## 16 regions with no links:
## 19 112 113 119 120 121 122 123 124 125 126 128 129 130 131 132
## 18 disjoint connected subgraphs
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9 10
## 16 3 4 12 19 19 25 16 11 6 1
## 3 least connected regions:
## 39 109 111 with 1 link
## 1 most connected region:
## 100 with 10 links
UTAM <- UTAM[!(UTAM$fid %in% c(19, 112, 113, 119, 120, 121, 122, 123, 124, 125, 126, 128, 129, 130, 131, 132)),]
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
#moran.test(UTAM$dist, listw = nb2listw(nb))
#UTAM$dist_leisure[is.na(UTAM$dist_leisure)] <- 0
moran <- c(0*(1:4))
for(i in 1:4){
m <- moran.test(data.frame(UTAM[,58+i])[,1], listw = nb2listw(nb))
moran[i] <- round(as.numeric(m$estimate[1]), digits = 2)
}
# Displaying the Moran's I
df <- data.frame(cbind(as.numeric(moran),motivo = c("Total", "Trabajo", "Estudio", "Esparcimiento")))
ggplot(df) +
theme_classic() +
geom_point(aes(x = motivo, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5) +
labs(x="Trip purpose", y="Moran") +
theme(legend.position = "left",
panel.border = element_rect(colour = "black", fill=NA, linewidth =0.5),
)
UTAMCentroids <- st_centroid(UTAM)
UTAMCentroids <- as(UTAMCentroids, "Spatial")
UTAMCentroids.geodata <- as.geodata(UTAMCentroids, data.col = "dist")
vario.ex<- variog(UTAMCentroids.geodata, bin.cloud=TRUE, option = "bin")
## variog: computing omnidirectional variogram
plot(vario.ex, main = "Semivariogram of the computed daily distance", cex.main = 1)
lines(vario.ex, type ="l", lty = 2, col="red", add = TRUE)
UTAM$Bogota <- 1
contour <- UTAM %>%
group_by(Bogota) %>%
summarise(geometry = st_union(geometry))
#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,100000,200000,500000,1000000,2500000,5000000)
# Defining the edge of Bogotá as the extent of the map
Emprise = as.owin(st_buffer(UTAM,0))
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 100000","100000 - 200000","200000 - 500000","500000 - 1000000","1000000 - 2500000","> 2500000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Total daily PKT per \nUTAM of residence", fill = "Total passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,20000,50000,100000,200000,500000,2000000)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_work)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 20,000","20,000 - 50,000","50,000 - 100,000","100,000 - 200,000","200,000 - 500,000","> 500,000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Commuting to work", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,10000,20000,50000,100000,250000,500000)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_school)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 10000","10000 - 20000","20000 - 50000","50000 - 100000","100000 - 250000","> 250000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Going to school or university", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
colors <- c("#f9f5a7", "#f7eca0", "#f3d893", "#e9b179", "#d9805a","#c64c3e")
bks = c(0,5000,10000,20000,30000,60000,100000)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_leisure)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 5000","5000 - 10000","10000 - 20000","20000 - 30000","30000 - 60000","> 60000"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Leisure or physical activity", fill = "Passengers-km per day", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM$dist_walk[is.na(UTAM$dist_walk)] <- 0
#Creating a color palette
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,0.5,0.8,1.2,1.5,2,4)
bks <- c(bks_walk,100)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_walk/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 0.63","0.63 - 0.75","0.75 - 1.03","1.03 - 1.32","1.32 - 1.58","> 1.58"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Walking", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM$dist_auto[is.na(UTAM$dist_auto)] <- 0
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,1,2,3,4,5,10)
bks <- c(bks_car,100)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_auto/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 0.56","0.56 - 1.01","1.01 - 1.86","1.86 - 3.02","3.02 - 4.65","> 4.65"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Private car", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM$dist_informal[is.na(UTAM$dist_informal)] <- 0
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
bks = c(0,0.2,0.5,1,2,3,10)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_informal/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 0.2","0.2 - 0.5","0.5 - 1","1 - 2","2 - 3","> 3"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Informal transport", fill = "Average daily distance per capita (km)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM$dist_BRT[is.na(UTAM$dist_BRT)] <- 0
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
bks = c(0,1,2,4,6,8,15)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_BRT/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_sf(data = Transmi, aes(color = "BRT (Transmilenio)"))+
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 1","1 - 2","2 - 4","4 - 6","6 - 8","> 8"))+
scale_color_manual(values = "red")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Bus Rapid Transit (BRT)", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
UTAM$dist_transit[is.na(UTAM$dist_transit)] <- 0
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "#bb2b30"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#efeff3", "#dddee9", "#caccde", "#aeb1ce", "#6c77a9","#284689")
#bks = c(0,2,4,6,8,10,20)
bks <- c(bks_transit,100)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$dist_transit/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 3.09","3.09 - 3.97","3.97 - 5.87","5.87 - 7.60","7.60 - 9.24", "> 9.24"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "Public transport", fill = "Average daily distance per capita (km)", linetype = "", color = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
viajes_dist <- viajes_dist %>% left_join(EF, by = "modo_principal")
viajes_dist$`CO2-eq` <- viajes_dist$`CO2-eq`*viajes_dist$dist_real/1000
viajes_dist$CO <- viajes_dist$CO*viajes_dist$dist_real/1000
viajes_dist$NOx <- viajes_dist$NOx*viajes_dist$dist_real/1000
viajes_dist$SO2 <- viajes_dist$SO2*viajes_dist$dist_real/1000
viajes_dist$COV <- viajes_dist$COV*viajes_dist$dist_real/1000
viajes_dist$PM.2.5 <- viajes_dist$PM.2.5*viajes_dist$dist_real/1000
viajes_dist$PM.10 <- viajes_dist$PM.10*viajes_dist$dist_real/1000
emissions_mode <- viajes_dist[,c(9,10,18,21:27)] %>%
group_by(modo_principal) %>%
summarize(`CO2-eq` = sum(f_exp*`CO2-eq`), CO = sum(f_exp*CO), NOx = sum(f_exp*NOx), SO2 = sum(f_exp*SO2), COV = sum (f_exp*COV), PM.2.5 = sum(f_exp*PM.2.5), PM.10 = sum(f_exp*PM.10))
#CO2 emissions are converted in metric tons.
emissions_mode$`CO2-eq` <- emissions_mode$`CO2-eq`/1000000
#CO emissions are converted in metric tons.
emissions_mode$CO <- emissions_mode$CO/1000000
#NOx emissions are converted in metric tons.
emissions_mode$NOx <- emissions_mode$NOx/1000000
#PM 2.5 emissions are converted in kg.
emissions_mode$PM.2.5 <- emissions_mode$PM.2.5/1000
#CO2 emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -`CO2-eq`), `CO2-eq`))+
theme_bw()+
geom_col(fill = 'darkorange2')+
labs(x="Main mode", y = "Greenhouse gas emissions (tCO2-eq/day)")+
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))
#CO emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -CO), CO))+
theme_bw()+
geom_col(fill = 'darkgreen')+
labs(x="Main mode", y = "Carbon monoxyde emissions (tCO/day)")+
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))
#NOx emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -NOx), NOx))+
theme_bw()+
geom_col(fill = 'royalblue')+
labs(x="Main mode", y = "Nitrogen oxydes emissions (t/day)")+
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))
#PM 2.5 emissions
ggplot(emissions_mode, aes(reorder(modo_principal, -PM.2.5), PM.2.5))+
theme_bw()+
geom_col(fill = 'darkred')+
labs(x="Main mode", y = "PM 2.5 Emissions (kg/day)")+
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=1))
UTAM <- st_read(dsn = "Data", layer = "EMU2019_area") %>% st_transform(4326) # UTAM boundaries
## Reading layer `EMU2019_area' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N
UTAM$Bogota <- 1
contour <- UTAM %>%
group_by(Bogota) %>%
summarise(geometry = st_union(geometry))
Emissions_UTAM <- viajes_dist[,c(1,9,21)] %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(CO2 = sum(f_exp*`CO2-eq`/1000000))
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")
bks = c(20,50,70,100,150)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = CO2))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkorange2")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "GHG Emissions per UTAM of residence", fill = "Total emissions (tCO2-eq/day)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Emissions_UTAM <- viajes_dist[,c(1,9,22)] %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(CO = sum(f_exp*CO/1000000))
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")
bks = c(1,2,3,5,10)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = CO))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkgreen")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "CO Emissions per UTAM of residence", fill = "Total emissions (tCO/day)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Emissions_UTAM <- viajes_dist[,c(1,9,23)] %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(NOx = sum(f_exp*NOx/1000000))
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")
bks = c(0.2,0.5,1,1.5,2)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = NOx))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#ebf7ff", high = "royalblue")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "NOx Emissions per UTAM of residence", fill = "Total emissions (tNOx/day)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
Emissions_UTAM <- viajes_dist[,c(1,9,26)] %>%
left_join(Manzana[,c(3,9)], by = c("id_hogar" = "ID_HOGAR")) %>%
group_by(UTAM)%>%
summarize(PM.2.5 = sum(f_exp*PM.2.5/1000))
UTAM <- UTAM %>% left_join(Emissions_UTAM, by = "UTAM")
bks = c(1,2,5,10,18)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = PM.2.5))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkred")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "PM.2.5 Emissions per UTAM of residence", fill = "Total emissions (kg PM 2.5/day)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#bks = c(0.5,0.6,0.8,1,1.2)
bks_co2 <- round(as.numeric(quantile(1000*UTAM$CO2/UTAM$NUMPERSTOT, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- bks_co2
UTAM$CO2[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = 1000*CO2/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkorange2", na.value = "grey80")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "GHG Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (kg CO2-eq)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
bks = c(30,40,50,60,80)
UTAM$CO[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = 1000000*CO/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkgreen", na.value = "grey80")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "CO Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (g CO)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#bks = c(1,2,3,4,6)
UTAM$NOx[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = 1000000*NOx/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(low = "#ebf7ff", high = "royalblue", na.value = "grey80")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "NOx Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (g NOx)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
bks = c(60,80,90,120,150)
UTAM$PM.2.5[UTAM$UTAM %in% c("UTAM3", "UTAM2")] <- NA
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = 1000000*PM.2.5/NUMPERSTOT))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_steps(breaks = bks, low = "#fafaaa", high = "darkred", na.value = "grey80")+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "PM.2.5 Emissions per capita \nper UTAM of residence", fill = "Average daily emissions per capita (mg PM 2.5)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
summary(nb)
## Neighbour list object:
## Number of regions: 132
## Number of nonzero links: 610
## Percentage nonzero weights: 3.500918
## Average number of links: 4.621212
## 16 regions with no links:
## 19 112 113 119 120 121 122 123 124 125 126 128 129 130 131 132
## 18 disjoint connected subgraphs
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9 10
## 16 3 4 17 19 18 29 8 11 6 1
## 3 least connected regions:
## 39 109 111 with 1 link
## 1 most connected region:
## 100 with 10 links
UTAM <- UTAM[!(UTAM$fid %in% c(19, 112, 113, 119, 120, 121, 122, 123, 124, 125, 126, 128, 129, 130, 131, 132)),]
nb <- poly2nb(pl = UTAM, snap = 50, queen = TRUE)
#moran.test(UTAM$dist, listw = nb2listw(nb))
contour <- UTAM %>%
group_by(Bogota) %>%
summarise(geometry = st_union(geometry))
UTAM$CO2[is.na(UTAM$CO2)] <- 0
UTAM$CO[is.na(UTAM$CO)] <- 0
UTAM$NOx[is.na(UTAM$NOx)] <- 0
UTAM$PM.2.5[is.na(UTAM$PM.2.5)] <- 0
moran <- c(0*(1:4))
for(i in 1:4){
m <- moran.test(data.frame(UTAM[,58+i]/UTAM$NUMPERSTOT)[,1], listw = nb2listw(nb))
moran[i] <- round(as.numeric(m$estimate[1]), digits = 2)
}
# Displaying the Moran's I
df <- data.frame(cbind(as.numeric(moran),motivo = c("GEI", "CO", "NOx", "PM 2.5")))
ggplot(df) +
theme_classic() +
geom_point(aes(x = motivo, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5) +
labs(x="Emission type", y="Moran") +
theme(legend.position = "left",
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
)
#UTAM must now be in projected coordinates
UTAM <- UTAM %>% st_transform(32618)
#Creating a color palette
colors <- c("#FAFAAA", "#F7DF88", "#F5C565", "#F2AA43", "#F09021", "#EE7600")
#bks = c(0,0.5,0.6,0.8,1,1.2,1.5)
bks <- c(bks_co2,50)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000*UTAM$CO2/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 0.48","0.48 - 0.55","0.55 - 0.69","0.69 - 0.87","0.87 - 1.07","> 1.07"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "GHG Emissions per capita", fill = "Average daily emissions per capita (kg CO2-eq)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "darkgreen"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#FAFAAA", "#C8DB88", "#95BE65", "#63A043", "#318221", "#006400")
bks = c(0,30,40,50,60,80,100)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$CO/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 30","30 - 40","40 - 50","50 - 60","60 - 80","> 80"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "CO Emissions per capita", fill = "Average daily emissions per capita (g CO)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
#pal <- colorRampPalette(colors = c("#ebf7ff", "royalblue"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#EBF7FF", "#C9DAF9", "#A6BEF3", "#84A1ED", "#6285E7", "#4169E1")
bks = c(0,4,6,8,10,12,20)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$NOx/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 4","4 - 6","6 - 8","8 - 10","10 - 12","> 12"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "NOx Emissions per capita", fill = "Average daily emissions per capita (g NOx)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))
#Creating a color palette
#pal <- colorRampPalette(colors = c("#fafaaa", "darkred"), interpolate = "linear")
#colors <- pal(6)
colors <- c("#FAFAAA", "#E3C888", "#CD9565", "#B76343", "#A13121", "#8B0000")
bks = c(0,60,80,90,120,150,180)
# Creating a ppp object (spatstat) and including the extent and values to be smoothed (i.e. the multiplication factor at 12.00 pm)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = 1000000*UTAM$PM.2.5/UTAM$NUMPERSTOT)
# Computing the smoothed surface (smoothing bandwidth: 1 km and picture spatial resolution: 1 ha) --> takes some time
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 800, weights = UTAM.ppp$marks, eps=sqrt(10000))
# Coercing the smoothed surface into a raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # to specify a CRS to the raster object
# Reclassifying the smoothed surface
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
# Vectorizing the reclassed surface
cartelissee.vecteur <- st_as_sf(st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# Sorting the layer column information
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ] %>% st_transform(4326)
# Configuring the window margins to maximize the extent of the map
par(mar = c(0, 0, 0, 0))
ggplot()+
theme_bw()+
geom_sf(data = cartelissee.vecteur, color = NA, aes(fill = factor(layer)))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_sf(data = contour, colour = "grey30", fill = NA) +
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 1.5, check_overlap = FALSE) +
scale_fill_manual(values = colors, labels = c("< 60","60 - 80","80 - 90","90 - 120","120 - 150","> 150"))+
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(tag = "PM 2.5 Emissions per capita", fill = "Average daily emissions per capita (mg PM 2.5)", linetype = "") +
theme(plot.tag = element_text(size = 9, face = "bold", hjust = 0), plot.tag.position = c(0.65, 0.9),
legend.title = element_text(size = 8),
legend.text = element_text(size=8),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(0.5, "cm"), width = unit(0.5, "cm"), style = north_arrow_orienteering(text_size = 4))